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ABSTRACT 

We present high resolution (R ~ 50,000) spectroastrometry over the CO P' overtone 
bandhead of a sample of seven intermediate/massive young stellar objects. These are primar- 
ily drawn from the red MSX source (RMS) survey, a systematic search for young massive 
stars which has returned a large, well selected sample of such objects. The mean luminosity 
of the sample is approximately 5x lO'* Lq, indicating the objects typically have a mass of ~ 15 
Mq. We fit the observed bandhead profiles with a model of a circumstellar disc, and find good 
agreement between the models and observations for all but one object. We compare the high 
angular precision (0.2 - O.SxlO"-' arc-seconds) spectroastrometric data to the spatial distri- 
bution of the emitting material in the best-fitting models. No spatial signatures of discs are 
detected, which is entirely consistent with the properties of the best-fitting models. Therefore, 
the observations suggest that the CO bandhead emission of massive young stellar objects orig- 
inates in small-scale disks, in agreement with previous work. This provides further evidence 
that massive stars form via disc accretion, as suggested by recent simulations. 



Key words: stars: early-type - stars: formation - stars: circumstellar matter 
accretion discs - line: profiles - techniques: high angular resolution 



accretion. 



^ 



1 INTRODUCTION 

Massive stars (My, > 8 Mq) can have a significant impact upon their 
environment, despite the fact they are less numerous than their 
lower mass counterparts. Massive stars can dominate their host 
galaxy's luminosity and inject prodigious amounts of energy into 
the interstellar medium (ISM). This injection of energy can regu- 
late subsequent star formation, and provides a key source of heat- 
ing and turbulence in the ISM. Furthermore, the enriched material 
injected into the ISM by supernovae accompanying the demise of 
massive stars forms a crucial component in subsequent generations 
of stars and planets. Theref ore, massive stars are imp ortant from 
galactic to planetary scales jZinnecker & YorkduOOTl) . However, 
despite their importance, our knowledge of how massive stars form 
is less complete than in the case of solar mass stars. 
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There has been considerable theoretical uncertainty over the 
formation of massive stars. A young massive star is expected to 
attain the luminosity of a main sequence OB star while it is still ac- 
creting material (according to a n extrapolation of the standard low 
mass star formation scenario of lShu et al.lll987h . As a result, it has 
been thought that the immense luminosity of massive stars could 
provi de sufficient radiation pressure to reve r se the in-fall of ma- 
terial JLarson & Starrfieldll97ll : lKahiilll974[Wolfire & Cassinellil 
Il987h . Consequently, alternative modes of massive star formation 
have been p roposed, for example competitive accretion a nd stel- 
lar mergers JBonnell et id]|l998l ; ISallv & Zinnecke3 12005^. How- 
ever, recent 3D hydrodynamic simulations demonstrate that radi- 
ation pressure does not prevent disc accretion forming stars of at 
least ~ 50 Mo JKrumholz et al.l2009l) . The key detail being that ac- 
cretion is confined to an equatorial disc, shielding the accreting 
material from the brunt of the radiati on, and channeling the ra- 
diation pressure into the polar regions jYorke & Sonnhaltenl2002l : 
iKrumholz et al.ll2009l : r\iidva et alj20m 

However, relating theoretical models to observations is chal- 
lenging. The Kelvin-Helmholtz timescale, the time taken for a 
proto-star to convert its potential energy to thermal energy and be- 
gin nuclear fusion, is approximately lO'* years for a massive star. 
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compared to ~ 10'' years for a star of 1 Mq. This short timescale, in 
conjunction with the innate rarity of massive stars, makes it difficult 
to catch a massive star in the act of forming. In addition, this short 
time scale means that massive stars join the main sequence (i.e. be- 
gin nuclear fusion in the core) before their natal material is cleared. 
Therefore, massive stars form under many magnitudes of extinc- 
tion. Finally, as massive stars are rare, they tend to be further away 
than nearby sites of low mass star formation. Consequently, study- 
ing massive young stellar objects (MYSOs) on scales of several au 
requires high angular resolution tec hniques such as spectroastrome- 
try an d optical interferometry (e.g. lOavies et al.ll20ld : lde Wit et al.l 
l2009l) . 

As a result, direct evidence for accretion discs around MYSOs 
- which we define as mid infrared bright stellar objects that have at- 
tained the luminosity of an early OB star but have yet to ionise their 
surroundings to form an Hn region - is sparse. Indirect evidence 
that massive stars form by disc accretion is provided by the dis- 
covery o f collimated jets and small scale outflo ws emanating from 
MYSOs JDavis et al.ll2004 IPavies et alJIlOlOl) . In a few isolated 
cases, direct detections of discs and flattened structures around 
MYSOs have been made ds hepherd ^TaLllOOll : iBeltran et al.l2005l : 
IPatel et al.l2005luOkamoto et al.i2009f) . However, such studies typ- 
ically probe dust emission at long wavelengths, and thus large dis- 
tances (of order hundreds to thousands of au) from the central star. 
Indeed, the rotating structures detected bv lBeltran et al.l ( I2005I) are 
too large to be considered circumstellar dis cs, as they are unsta - 
ble and may fragment (see the discussion in lCesaroni et aLluOOTT) . 
Therefore, to prove that MYSOs form via disc accretion requires 
evidence for small scale, au sized, gaseous discs. 

A powerful diagnostic of small-scale discs around young 
stellar objects (YSO)s is provided by CO 1^' overtone band- 
head emission at 2.3 yum. Such emission requires high densities 
(n > lO'^cm"^) and high temperatures (T~ 2 500K), ind icative of 
disc material close to the stellar surface (e.g. ICarjl 19891) . In addi- 
tion, such emission associated with YSOs can be well fit with a 
model of the band head emission originating in a Ke plerian disc, a 
few au in size (e. g. ICarnl 19891: Chandler et alJI 19951) . In particular, 
iBik & Thil J2004h and lBlum et alj ( l2004h found that the CO band- 
head profiles of intermediate and massive YSOs could be fit with 
emission from discs with small, » 0. 1 au, inner radii. This provides 
strong support for the formation of massive stars via disc accretion. 
However, the number of objects observed so far is small, and there 
are several intermediate mass objects among the observed samples. 

In this paper we extend the work of iBik & Thil j2004) and 
iBlum et all ( |2004|) to a sample of MYSOs drawn from the red 
MSX (Midcourse Space ex periment) source (RMS) survey cata- 
logue JUrguhart et al.ll2008h . The RMS is a galaxy wide survey of 
MYSOs and the most representative sample of this class of objects. 
Previous se arches for MYSOs relied on the IRAS Point Source Cat- 
alogu e (e.g. lChan et"ai]|l996l : lM"olinari et al.lll996l : ISridharan et al.l 
l2002h . However, due to the large beam of IRAS (2-5' at 100 ^m), 
such studies suffered from considerable source confusion, and 
were biased away from the Galactic Plane, where the majority of 
MYSOs are expected. The RMS, however, utilises the MSX sur- 
vey of the galactic plane in the mid infrared l lPrice et alj200lh . The 
MSX survey has a resolution of 20" and thus offers a factor of 50 
improvement in spatial resolution over IRAS , which allows sources 
to be detected in regions that are otherwise unresolved. Therefore, 
the use of the MSX survey allows a unique and representative 
sample of MYSOs to be selected. Specifically, the RMS survey 
used colour selection criteria and the MSX and 2MASS catalogues 



JEgan et alj |2003|; ICutri et alj 120030 to select an unbiased sam- 
ple o f approximately 2000 potential MYSOs (see iLumsden et al.l 
[2003). However, this initial sample contains objects such as plan- 
etary nebulae, Hn regions and low luminosity YSOs that have a 
similar appearance to MYSOs in the near and mid infrared. These 
contaminant objects have been eliminated via an extensive multi- 
wavelength campaign featu ring high resolution (1-2" ) observa- 
tions in the radio continu um JUrguhart et al.'2007'.'200^. the '^CO 
J=l-0 and J=2-l lines JUrguhart et al. 2007. 2008) the mid in- 
frared JMottram et al.l2007l) and the near infrared (e.g. lClarke et al.1 
I2OO6I) . In total, the RMS databastQ provides a large (-500), well- 
selected sample of mid infrared bright MYSOs. 

Here, we exploit the RMS survey to study the accretion char- 
acteristics of a sample of MYSOs. From the low resolution, NIR 
spectroscop y undertaken as a p art of the classification stage of the 
survey (e.g. Iciarke et alJuOOo) . we selected a sub-sample of ob- 
jects with CO P' overtone bandhead emission at 2.3 /jm. We have 
observed these objects at high spectral resolution to investigate 
whether their CO bandheads are consistent with emission originat- 
ing in a small-scale circumstellar disc. 

To fully constrain the models, we also require information on 
the spatial distribution of the emission. However, the emission re- 
gion is expected to be small. At a typical distance of Ikpc, a disc 
of 1 au in size subtends an angle of only 1 milli-arcsec (mas). 
Spectroastrometry is one of a few approaches that offers the re- 
quired, sub-mas, angular precision and high spectral resolution. 
Indeed, such an approach has already been shown to be able to 
detect and characteris e circumstellar discs wit h a precision of ap- 
proximately 0. 1 mas jPontoppidan et al.ll2008l . Wheelwright et al. 
in prep.). Therefore, we use spectroastrometry to probe the spatial 
behaviour of the bandhead emission with high spectral resolution 
(which is required to obtain kinematic information). 

The paper is structured as follows: in Section[2]we present the 
sample of MYSOs, details of the observations and the data reduc- 
tion processes used. The resultant data are presented in Section [3] 
alongside the results of fitting a model of CO emission from a cir- 
cumstellar disc to the data (the model is described in Section [T2] ). 
Section|4]presents a discussion of the results. We conclude the pa- 
per in Section|5] 



2 OBSERVATIONS AND DATA REDUCTION 

2.1 Observations 

The sample was selected using the low resolution A^-band spec- 
troscopy undertaken as part of the RMS survey (Cooper et al. in 
prep.) to identify objects with CO P' overtone bandhead emission. 
In addition to the resulting targets, we included IRAS 08576-4334 
and M8E in the sample. IRAS 08576-4334 is an intermediate 
mass YSO known to exhibit CO first overtone emission JBik & Thil 
I2OO4JD. This object was included in the sample of MYSOs as its 
CO e mission has been repor ted to be extended on scales of tens 
of au JGrave & KumajuOOTl) . making it an appropriate target to 
study the spatial distribution of CO emission. M8E is a very bright 
MYSO that is thought to possess a circumstellar disc JSimon et al.l 
|i985D. The final sample is presented in Table [Tj along with details 
of the observations. 
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High resolution spectroscopy of the sample at 2. Sum wa s ob- 
tained usi ng Phoenix at Gem ini South jHinkle et alj 120001) and 
CRIRES l lKaeufl et alJ |2004|) on UTl at the VLT. When using 
Phoenix a slit of 0.34 arcsec was used, which resulted in a spec- 
tral resolution of approximately 50 000, or 6 km s ' . The CRIRES 
observations were conducted using a slit of 0.4 arcsec, and the same 
spectral resolution. Observations were conducted using 4 slit posi- 
tion angles (PA): 0°, 90°, 180° a nd 270° , whi ch is a requirement for 
accurate spectroastrometry (see lBailevlll998h . 

The observations with Phoenix were conducted with the 
K4396 filter, which has a central wavelength of 2.295 //m, to isolate 
the region surrounding the CO P' overtone bandhead. When us- 
ing CRIRES, the Kj filter was used with the spectral configuration 
identified by the reference wavelength 2.2932 ^m. This resulted 
in the CO P' overtone bandhead being located on the third chip, 
with a few ro- vibrational lines evident on the fourth (e.g. between 
~ 2.30 - 2.31/jm). All observations were conducted using a stan- 
dard nodding sequence along the slit to remove the sky background. 
Where a natural guide star was available (G332. 8256-00. 5498 and 
G347.0775-00.3927) the adaptive optic capabilities of the VLT, the 
Multi-Application Curvature Adaptive Optics facility, was utilised. 
All the observations with Phoenix at Gemini South were con- 
ducted with natural seeing. The resulting Full Width at Half Maxi- 
mum (FWHM) of the spectral spatial profiles ranged from 0.27 to 
0.75 arcsec, and was typically 0.5 arcsec. 

Telluric standard stars, late B type or early A type dwarfs, 
were observed with the identical instrumental setup as the science 
observations, and at similar air-masses. These spectra were subse- 
quently used to connect the science spectra for telluric lines, and to 
provide a rough photometric calibration. 

2.2 Data reduction 

The data were reduced in a standard fashion. A master flat frame 
was constructed, connected for dark current and normalised. The in- 
dividual exposures were then corrected by the normalised, average 
flat frame (the IRAS 08576-4334 and G287.3716-f00.6444 data 
were corrected with a median smoothed flat frame due to a dis- 
crepancy between the flat response and the behaviour of the A and 
B spectra). Pairs of spectra defined by the A and B nodding posi- 
tions were subtracted from one-another, and the individual intensity 
spectra were extracted from the resultant data. Spectroastrometry 
was conducted by fitting a Gaussian profile to the spatial profile of 
the spectrum at each dispersion pixel. This resulted in a position 
spectrum; the photo-centre of the spectrum as a function of wave- 
length, associated with each conventional spectrum. The centroid 
may be determined with high precision, allowing the centre of the 
emission to be traced to within a fraction of a pixel (for example see 
Oudmaijer et al. 2008; van der Plas et al. 2009; Wheelwright e t al.l 
20101) . Visual inspections were used to ensure the Gaussian pro- 
file was a va lid representation o f the Point Spread Function. As 
noted by e.g. lTakami et alj l l2003h . the precision of spectroastrom- 
etry scales linearly with the width and signal to noise ratio (SNR) 
of the spatial profile. As a consequence of the high SNR data (typi- 
cally 200-300), and narrow spatial profiles, we were able to achieve 
an angular precision down to approximately 0.2 mas. 

The individual flux spectra at a given PA were combined, as 
were the positional spectra, to form an average spectrum for the 
PA in question. The intensity spectra were divided by the spec- 
trum of a telluric standard to correct for telluric absorption lines. 
Finally, the intensity spectra at each PA were combined to form 



the average for the object in question. Dispersion correction was 
performed using the telluric lines in the telluric standard spectra 
and the high resol ution NIR spectrum of Arcturus presented by 
iHinkle et al.l jl995r) . The dispersion correction typically had an rms 
of 0.5 km s"' , i.e. less than -j^ the resolution. The positional spectra 
at anti-parallel position angles were combined to form the average 
North-South and East- West position spectra as follows: (0 -180)/2 
and ( 90-270)/2. This eliminates artificial signatures (see iBailevi 
Il998l) . which do not follow the slit rotation. 

Finally, the object spectra were flux calibrated using the tel- 
luric standards. Dividing the science data by a telluric standard 
spectrum allowed the continuum brightness ratio to be determined. 
The A^-band magnitude of the science object was then be esti- 
mated from the flux ratio and the telluric's A'-band magnitude from 
2MASS. The result was generally within 10 per cent of the pub- 
lished A'-band magnitude. Therefore, the objects' A"-band magni- 
tudes were used to estimate the continuum flux density. The flux 
over the bandhead was then determined based on the observed 
strength of the CO emission. 



2.3 Extinction determination 

In Section [331 we compare the flux densities of the best-fitting mod- 
els with the observed values, a consistency check that is often not 
performed in similar studies. However, in order to assess the intrin- 
sic bandhead flux, we must first determine the extinction towards 
each object. 

To achieve this we u se the methodology oflPorter et alj l ll998h . 
the NIR extinction law of lStead & Hoard ( 120090 and the low resolu- 
tion, /f -band spectroscopy undertaken as part of the RMS survey. It 
is important to note that the general extinction law for the interstel- 
lar medium may not be applicable to regions of high extinction (e.g. 
see lMoore et alj|2005j). However, the fla ttening of the NIR extinc- 
tion law reported bv lMoore et al.l ( 120051) only becomes apparent at 
an extinction of approximately 4 magnitudes in K, an d in general 
the ex tinction towards MYSOs is less than this (e.g. IPorter et al.l 
Il998l) . Therefore, the general interstellar extinction law should be 
applicable to the sample presented here. 

It is assumed that the science targets are relatively hot, early 
type stars. The NIR continuum of a hot star is we ll approximated 
by the Rayleigh-Ieans (RJ) tail jPorter et alJI 19981) . By reddening 
the RJ slope, Fj oc A^*, until it fits the observed //-band contin- 
uum, the extinction can be determined. This is shown in Fig.[T]for 
the RMS sources for which we have low resolution NIR spectra. 
Only the //-band spectral continuum is used, as it is less likely 
to be contaminated with emission from hot dust than the A'-band. 
Extinction estimates for MYSOs derived from the continuum in 
the //-band are generally consistent with the results of other meth- 
ods, s uch as using hydrogen recombination line ratios IPorter et al.l 
Il998r) . However, extinction values based on the continuum from 
the //-band and into the A'-band tend to be overesti mates, due to 
infrar ed excess emission from hot circumstellar dust dPorter et al.l 
119981) . Therefore, we use the wavelength range ~ 1.5 - 1.6 yum to 
estimate the extinction towards the science targets. The determined 
values of Av and Af^ are presented in Table|2] 

The uncertainty in determining the extinction towards an ob- 
ject u sing the slope of th e its NIR continuum is typically 10-15 per 
cent dPorter et al.ll 19981 based on visually assessing the quality of 
the fit). We find the uncertainty based on the x^ statistic associated 
with fitting the continuum is also generally 10 percent, or approxi- 
mately 0.1-0.3 magnitudes. While this uncertainty is in agreement 
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T able 1. A summary of the sample and observations. The bolometric luminosities (Col. 5) are taken from the RMS database' , which are based on the work 
of iMottram et alj 1201(1) and Mottram et al. (2010: in prep.). The kinematic distances (Col. 7) are also taken from the RMS d atabase", unless tak en from 
the literature as indicated. The masses (Col. 4) are estimated from the bolometric luminosities and main sequence relationships jMartins et alj |2005j). unless 
otherwise stated. The X'-band magnitudes (Col. 6) are from the 2MASS Catalog^ , unless otherwise stated. The integration time (Col. 8) denotes the total 
integration time per object, and the precision of the re-binned spectroastrometric data (the average rms of the position spectra) are presented in Col. 9. Finally, 
we list the dates the objects were observed (Col. 10). 



Name 



RA 

J2000 



Dec 
J2000 



M 

Mo 



LboI 



K 

mag 



kpc 



Time Apos Date 
hours mas 



Gemini South sample: 

IRAS 08576-4334* 

G287.3716+00.6444 

M8E* 

VLT sample: 

G308.9176+00.1231 

G332.8256-00.5498 

G347.0775-00.3927 

G033.3891+00.1989 



08 59 25.2 
10:48:04.6 
18 04 53.3 

13:43:01.6 
16:20:11.1 
17:12:25.8 
18:51:33.8 



-43 45 46.0 
-58:27:01.0 
-24 26 42.3 

-62:08:51.3 
-50:53:16.2 
-39:55:19.9 
+00:29:51.0 



6.r 

15.4 
I3.5A 

42.6 
24.5 
17.7 
11.8 



9.4* 
2.9X10-* 7.5 

4.4 

3.9x10^ 6.4 

1.3x10^ 8.9 

4.3x10^* 8.5 

1.3x10'' 7.2 



2.2* 4.27 0.76 2008: 01/02, 01/10, 01/23 

5.6 3.20 0.25 2008:01/10,01/21 

1.9* 0.27 0.47 2007: 10/04 

5.3 0.60 0.13 2009:04/08,04/11 

3.8° 2.40 0.32 2009: 05/08, 08/19, 09/01, /09/08 

14.9 2.80 0.15 2009: 07/17, 07/19, /08/03, /08/15 

5.5 1.22 0.63 2009:05/14,07/05,08/02,08/15.08/19 



t: |http://www.ast.leeds.ac.uk/RMS/| o: the distai ices in the RMS datab ase are from lSronfman et all fT993):IUrauhart et al <2007ll2008l). t : Cutri et alj J2003h. 
• : Kin ema tic distance from the rotation curve of lBrand & Blitj h993h and a Vur o f 7.5km s~ ' JBronfman et alj 19961) , »: IChini & Neckell il98ll), *:lBiketalJ 
i200d) . /\: lLinzetai] i2009l) . o: based on the photometry andKvsJ-K diagram of lBik et alj ilOOdt) . a distance of 2. 17 kpc and the data of lHarmanecUl988h . 
n: Near/Far ambiguity, here, we use the smallest possible distance. 



with previous work there is an important caveat to consider. This 
approach assumes the continuum observed is purely photospheric 
in origin. If the stellar continuum is contaminated by hot dust emis- 
sion, the blue to red slope will become steeper, and the extinction 
will be overestimated. Conversely, if continuum is contaminated by 
scattered light the extinction will be underestimated due to a shal- 
lower slope. 

In general the values of Ay are similar, with most values 
close to 20-30. However, the derived extinction of Av=10 to 
G332.8256-00.5498 is significantly higher than the other values. 
Contamination of the continuum by dust emission may have re- 
sulted in the extinction being over estimated as the wavelength 
range used to determine the extinction was extended to 1.6-1.7 yum 
(as no flux was observed at 1.5-1.6 yum). Assessing the presence 
of dust-excess emission by SED modelling is beyond the scope 
of this paper. However, we can use the ratio of the Hi lines in 
the low resolution spectru m to assess the extinction independently 
of the continuum (see e. g. iLandini et al.ll 19841 : iLumsden & PuxlevI 
ll996l : lMoore et al.ll2005l) . Using the ratio of the By and BrlO lines 
and a value of 2.1 for the slope of the NIR extinction law re- 
sults in a value of Ak=2.4±0.3. This value is consistent with the 
mean value of the sample, and since it is not affected by emis- 
sion from hot dust is likely to be closer to the correct value than 
the previous value of 7.9. Therefore, we use this value of Ak for 
this object. We note that determining the extinction towards an 
object via the ratio of Hi recombination lines is only valid if the 
intrinsic lin e ratios can be estimated via hydrogen recombination 
models (e.g. Storey & Hummeill 19951) . This requires that case B of 
iBaker & Menzell ( 119381) applies, and in turn limits this method to 
Hii regions. As G332.8256-00.5498 is the only Hii region in the 
sample this is the only object we can apply this method to. 

These extinction estimates are subject to some uncertainty (of 
the order 10 per cent), which translates to a substantial uncertainty 
in the de-reddened continuum fluxes (up to approximately 50 per 
cent). None the less, the derived fluxes will provide a useful check 
on the modelling results. 



3 RESULTS 

3.1 The spectra 

Continuum normalised spectra of the MYSO sample at 2.3yum 
are presented in Fig. [2] The objects exhibit a range of CO band- 
head profiles. Some display a clear blue 'shoulder' at the band- 
head (e.g. G332.8256-00.5498), while others exhibit a relatively 
sharp rise from the blue continuum to the bandhead (e.g. IRAS 
08576-4334). The peak fluxes for the RMS objects are of order 10 
per cent the continuum, much less than that of IRAS 08576-4334 
(~ 80percent). Bandhead profiles with a prominent b lue shoulder 
are in dicative of emission from a rotating disc (e.g. iNaiita et al.l 
1 19961) . None the less, the profiles with sha rper edges to the band- 
head can also be fit with disc models (see iBik & Thill2004l) . pro- 
vided the disc is large or the inclination is low, thereby minimising 
the rotational broadening of the bandhead. 

We used a simple model to first fit the bandhead profiles. We 
then compared the predicted spectroastrometric signatures of the 
best-fitting models with the data. 



3.2 The model 

The model of CO P' overtone bandhead emission is based on 
iKraus et al.l OOOOh . Initially, a simple, geometrically flat disc is 
constructed, within which the excitation temperature and surface 
number density decrease with increasing radius. The decrease in 
temperature and surface number density are treated as power laws. 
The temperature decreases with R^"-''^ and the surface number den- 
sity falls off w ith /f"'^, in line with standa r d accr etion disc theory 
jPringlelflgsTh . and following [Kraus et all ( |2000|) . who also mod- 
elled CO emission from a circumstellar disc. As the disc model is 
very simple and does not incorporate accretion the temperature of 
the disc is determined as follows: T{R) = Teff x(^)"'"^, where T^ff 
is set based on the luminosity of the star (assuming it to be on the 
main sequence). The disc is split into radial and azimuthal cells. If 
the temperature in a cell is greater than 5000 K (the assumed de- 
struction temperature of CO), the flux from the cell is set to zero. 
If the outer radius of the disc is such that the outer cell tempera- 
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Figure 1. The RMS sources of which we have obtained low resolution NIR spectra as part of the follow-up program of observations for the RMS survey. The 
dashed line represents the Rayleigh- Jeans tail of a hot star (F,i oc i"'*). The solid line is this tail reddened to fit the slope of the //-band spectra. 

Table 2. The extinction (Cols 2 & 3) towards each object (Col. 1). The Ay values are determined from the best-fitting Ak values and the relationship between 
optical and NIR extinction from lCardelli et alj (j_989). In addition, we present the /f-band continuum flux density of each object, based on its K-brnid magnitude 
(Col. 4), and the estimated flux density at the peak of the CO P' overtone bandhead (Col. 5). Finally, we de-redden the CO bandhead flux (Col. 6) and compare 
it to the flux of the best-fitting model (Col. 7). The uncertainty in the extinction values are estimated from assessing the fit to the continuum with x^ (where 
the fitting range was specified to avoid lines and artefacts in the spectra). The uncertainty in the extinction is then used to determine the error in the flux. 
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tures are less than 1000 K the disc is shrunk until the temperature 
at the edge is 1000 K. Therefore, while the radii were varied during 
the fitting process, they are tied to the temperature structure of the 
disc. The CO e mission of each cel l is calculated according to the 
methodology of iKraus et alj OOOCi) . which is briefly described in 
the following. 



The population of the CO rotational levels (up to J=100) for 
the 2-0 vibrational transition in each cell is determined assuming 
local thermodynamic equilibrium. Then, the absorption coefficient 
is determined from the populations of the possible energy levels, 
and the transition probabilities of the respective lines. Assuming the 
absorption coefficient is constant along the line of sight, the optical 
depth is given by the product of the absorption coefficient per CO 
molecule and the CO column density. The column density is given 
by the surface number density, as we use a thin disc. The Dunham 
coe fficients re quired to calculate the CO energy levels were taken 
from lFarrenq e t al. ( 1991), and the Ein stein coefficients of ea ch ro- 
vibrational transition were taken from IChandra et alj ( Il99q) . The 
intrinsic line profile is assumed to be Gaussian, with the line-width 
being a free parameter. Once the individual cell spectra are calcu- 
lated, they are summed to create the total spectrum. This is then 



convolved with a Gaussian profile to match the spectral resolution 
of the observations. 

To determine the best-fitting model, the downhill simplex al- 
gorithm was used. This algorithm is supplied with the Interactive 
Data Language distribution (IDL), and is based on the amoeba rou- 
tine. The free parameters are: the inner and outer radius of the disc, 
the inclination of the disc, the surface number density at the in- 
ner edge of the disc and the line width within the disc. Once the 
best-fitting parameters have been determined the spectroastromet- 
ric signature of the best-fitting model is predicted. This is done by 
mapping the flux from each cell onto an array with orthogonal spa- 
tial and dispersion axes that represents a long-slit spectrum. The 
array is then convolved with a Gaussian profile in the spatial direc- 
tion, to represent the seeing conditions. The synthetic signature is 
then generated by fitting a Gaussian profile to the spatial distribu- 
tion at each dispersion pixel. 

The stellar mass is not a free parameter and is set to that listed 
in Table [T] which is generally determined fr om the luminosity o f 
the source and main sequence relationships JMartins et alJuOOSi) . 
However, the luminosities of IRAS 08576-4334 and M8E are not 
known. Therefore, the mass of IRAS 08576- 4334 was determ ined 
from its position in the Kys, ] -K diagram of lBik et al.l ( 120061) . We 



6 H. E. Wheelwright et al. 




2.293 2.294 2.295 2.296 2.297 
X (/im) 




2.293 2.294 2.295 2.296 2.297 
X (|iim) 



2.293 2.294 2.295 2.296 2.297 
A (/.im) 





1.10 


G30B.9176+00.1231 - 


■d 




L 


UJ 






a: 




/ ^rl.. 














(0 

e 




1.05 


if ' 


'"""^ 






J 


fc 


1.00 


f ; 



2.293 2.294 2.295 2.296 2.297 
A (/im) 



2.293 2.294 2.295 2.296 2.297 
A (/im) 




2.293 2.294 2.295 2.296 2.297 
A (/im) 



1.10 


j^G347.0775-00.39Z7 - 


E 


/ 


n 


Toil. . 





,/ 




ri 


fM 


t^ 1.00 


r - 


r 



2.293 2.294 2.295 2.296 2.297 
A (/im) 

Figure 2. The spectra of the sample around the CO P' overtone bandhead. The spectra have been divided by standard spectra, normalised, and re-binned 
by a factor of 3. The dashed line represents the best-fitting model, which is described in Section 13.21 The absolution features present in the spectra of 
G033.3891+00.1989 and G332.8256-00.5498 are residual telluric effects. 



revise the distance to IRAS 08576-4334, and find it is still located 
in the mid- to early B area of the K \& J - K diagram, and hence 
assume it h as the mas s of a B3 type star. The mass of M8E was 
taken from iLinz et alj | |2009|) . an d is the mass of the b est-fitting 
model from the grid of models bv lRobitaille et alj d2007h . 

Using the luminosity of a MYSO to determine its mass is sub- 
ject to some uncertainty as the observed luminosity may be due in 
part to accretion, resulting in an overestimate of the stellar mass. To 
evaluate the possible effect this may have, w e consulted the models 
of acc reting, massive protostars presented bv lHosokawa & Omukail 
( 120091) . During the early adiabatic accretion phase, the accretion 
luminosity of a massive protostar is greater than its intrinsic lu- 
minosity. However, by the time a massive protostar enters its main 
sequence accretion phase the intrinsic luminosity is greater than the 
accretion luminosity. While the resultant total luminosity is slightly 
greater than the zero age main sequence luminosity the two are not 



different by an order of magnitude. Given the high luminosity of 
the sample objects, it is likely they are in the Kelvin-Helmholtz 
contraction or main sequence accretion phases. Therefore, it is un- 
likely we significantly overestimate their luminosity. Furthermore, 
given that the objects are most likely in their main sequence ac- 
cretion phase, the main sequence relationship between mass and 
luminosity should be applicable. Therefore, it is surmised that us- 
ing main sequence relationships is sufficient for the purposes of this 
paper. 



3.3 Model fits and results 

Over the wavelength range at our disposal few individual ro- 
vibrational lines are visible. The fitting process is therefore dom- 
inated by the shape of the CO bandhead peak and shoulder. The 
appearance of the shoulder is set by the rotational profile of the in- 
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Figure 3. The spectroastrometric signatures of the sample over the CO l"' overtone bandhead emission, shown normalised to the continuum position. The 
data, the solid line, have been re-binned with a bin size of 65 km s"' (approximately 1 1 times the resolution element). We also present the predicted signatures 
of the best-fitting model (see Section [33t for each object, bar G287. 3716+00. 6444 (see Section l4!3t . to illustrate the signature of a rotating disc. The model 
signature is shown twice (at the un-binned spectral resolution): 1 . with its predicted amplitude (dashed line), and 2. with an amplitude of 2 mas in the negative 
direction to enhance its visibility (offset dot-dashed line). The maximum predicted signature is presented in the lower right comer of each plot (in mas). 



dividual lines, and is thus dependent on the inclination and the inner 
and outer disc radii (as the stellar mass is not varied). However, the 
temperature of the disc is also a function of distance from the cen- 
tral star. In general few individual lines are evident, indicative of 
relatively high temperatures - and/or high rotational broadening. 
As a result, the inner disc radii are generally small, of the order 
0. 1 au. Changing the outer radii has little effect on the quality of 
the fit as the surface number density, and thus flux, fall off steeply. 
As a consequence, it is the inclination that dominates the fit to the 
bandhead profile. The object with the most prominent blue shoul- 
der is fit with the largest inchnation (G332.8256-00.5498: ~ 42°), 
while the objects with the steepest slopes from the blue continuum 
to the bandhead peak (IRAS 08576-4334 and M8E) are fit with the 
lowest inclinations (~ 18° & 16°). 



The best fit parameters are presented in Table [3] and the re- 
sultant bandhead profiles are plotted with the data in Fig. [2] As 
can be seen from the reduced x~ values (the mean is 1.68) the 
model of CO emission originating from a circumstellar disc gen- 
erally provides a good fit to the data. The sole exception to this is 
G287. 3716+00.6444, which is discussed in more detail in Section 
|43] 

We compare the predicted bandhead fluxes with the un- 
reddened observed fluxes as a consistency check. As mentioned in 
Section [23l these un-reddened fluxes are inevitably subject to con- 
siderable uncertainties, yet still serve as an important validity test. 
We list both model and source fluxes in Table |2] These are similar 
to within a factor of 2-3, and generally consistent within the uncer- 
tainty from the extinction estimate (neglecting the addition error 
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due to uncertainties in the kinematic distances). Therefore, the ob- 
served and model fluxes are essentially consistent. This provides a 
further confirmation that disc models provide a good fit to the data. 



tra eliminates artefacts (the data of I Grave & Kuman ( 120071) were 
obtained using a constant PA.) 



3.4 Spectroastrometric signatures 

The spectroastrometric signatures associated with the spectra pre- 
sented in Fig.[2]are displayed in Fig.[3] The spectroastrometric data 
have been re-binned by a factor of ~ 11 times the resolution ele- 
ment. This rebinning factor was chosen as it was the maximum 
bin size that could still resolve the expected spectroastrometric sig- 
natures across the bandhead. The resulting average positional pre- 
cision is approximately 0.4 mas. In principle, this means we are 
probing au size scales at the kpc distances of the sample. Here, 
we compare the spectroastrometric data to the best-fitting models 
to determine whether these data can be used to further probe the 
circumstellar environments of the sample. 

The spectroastrometric signatures associated with the best- 
fitting models are also presented in Fig. [3] The predicted signa- 
tures are generally small, and to illustrate their appearance, an en- 
hanced spectrum is shown offset in each panel. The largest excur- 
sion occurs on the blue side of the bandhead, which is due to the 
asymmetry of the bandhead profile. Over the blue shoulder of the 
bandhead all the emission is blue-shifted. Therefore, at these wave- 
lengths the side of the disc coming towards the observer is signif- 
icantly brighter than the other side, which results in the prominent 
spectroastrometric signature at these wavelengths. Conversely, red- 
wards of the bandhead peak the receding side of the disc is bright- 
est, resulting in a positional excursion in the opposite direction. 
However, at these wavelengths red-shifted emission from the peak 
of the bandhead is mixed with blue-shifted emission of ro- vibration 
lines at longer wavelengths. As a consequence, the spectroastro- 
metric signature is smaller to the red of the bandhead peak than on 
the blue side. Over individual, resolved ro- vibrational lines the two 
sides of the disc are almost equally bright, and thus the excursions 
at positive and negative velocities are nearly s ymmetrical, result- 
ing in the classic 'S' shaped si gnature of a disc jPontoppidan et al.l 
l2008l : lvan der Plas et alj|2009l. Wheelwright et al. in prep.). As the 
individual ro- vibrational lines are located close together the adja- 
cent positional signatures result in an oscillating signature. 

The largest predicted excursions are of the order 0.01 mas, 
which is much smaller than the precision of the spectroastrometric 
data. The small scale of the spectroastrometric signatures is largely 
due to the large distances to the science objects, coupled with rel- 
atively weak CO emission, which is typically 10 per cent of the 
continuum flux. The largest predicted signatures are those of IRAS 
08576-4334 and M8E (both 0.017 mas), due to their proximity. 
Even with a positional precision of 0.15 mas, the best in the sam- 
ple, such excursions would not be detected. 

While the spectroastrometric data do not provide useful con- 
straints to the model fitting process, they are entirely consistent 
with the best-fitting models of emission originating in circumstel- 
lar discs. Moreover, the data do provide additional constraints on 
the source of the emission. We note that the spectroastrometric 
signature of IRAS 08576-4334 exhibits no detectable excursion 
with a pre cision of approximately 0.7 mas. This is contrary to the 
finding of IGrave & K umar ( 2007), who report an excursion of ap- 
proximately 18 mas, corresponding to a minimum size of 40 au 
at 2.17 kpc, and likely several times this. We suggest that the fea- 
ture in their data is an instrumental artefact, which is not present in 
our data as our practice of subtracting anti-parallel position spec- 



4 DISCUSSION 

4.1 A discussion on tlie best-fitting parameters 

The best-fitting inner disc radii are typically of the order 0.1 au, and 
the best-fitting outer disc r adii are generally 1-5 au. These values 
are consistent with those of lSik & Thil ( 120040 . and indicate that the 
CO emission arises from small scale circumstellar discs. At such 
small distances from the central star, CO molecules should be dis- 
associated by ultraviolet photons. However, the best-fitting number 
surface dens ities are sufficient for the CO molecules to self-shield, 
as found bv iBik & Thil J2004l) . We note that the simplistic treat- 
ment of the temperature gradient within the disc allows CO to be 
excited within 1 au of the central star. However, a full treatment 
of viscous heating in such a disc keeps t he mid-plane tempe rature 
above ~5000 K out to approximately 1 au dVaidva et al.l2009l ). This 
would have the effect of moving the inner radius of the CO emit- 
ting region further from the central star. The difference between 
the best-fitting inner radii and the radius at which the mid-plane 
temperature of a massive accretion disc drops below 5000K is of 
the order a few, and thus is unlikely to significantly affect the best- 
fitting properties. In addition, the viscous heating is dependent upon 
the accretion rate, which is uncertain and thus incorporating vis- 
cous heating in the model would introduce an additional unknown. 
Therefore, we judge the simplistic temperature gradient to be suffi- 
cient for the purposes of this paper. 

In the model the excitation temperature in the disc varies as a 
function of radius. If the disc is sufficiently large the outer tempera- 
ture falls below 1000 K and the outer cells contribute no flux to the 
total spectrum. In this situation increasing the outer radius will not 
affect the resultant ;^'^, as the output spectrum will not be changed. 
However, the physical disc may well be much larger than the CO 
emitting region, which constitutes the warmer and denser part of 
the disc. In this respect the outer r egions of the C O emitting region 
are consistent with the models of IVaidva et al.l ( I2009l) , who show 
that accretion discs around massive protostars may well be stable 
out to approximately 80 au. The dust sublimat ion radii associated 
with MYSOs are typically greater than 10 au dde Wit et aLll2007l: 
IVaidva et alj|2009l) . Therefore, the CO emission of the best-fitting 
models originates from within the dust sublimation radius, and thus 
offers a unique tracer of gas prior to its accretion onto the central 
star. 

The inner radii set the initial excitation temperature, which is 
a function of the distance from the central star. As a result, chang- 
ing the inner radius effects both the temperature in the disc and the 
rotational broadening of the bandhead. To consistently determine 
the temperature structure of the disc, and thus more stringent con- 
straints, a detailed radiative transfer model of the disc is required. 
However, even if the gas disc were modelled in detail, the proper- 
ties of the central star are essentially unknown. For example the star 
may swell to several ti mes its main sequence size as a consequence 
of high accretion rates dHosokawa & Omukail 20090 . which will re- 
duce its effective temperature. Therefore, even a sophisticated disc 
model wold be subject to several unknowns. 

The best-fitting line-widths are generally greater than the 
value of ~ 4kms"' determined bv lBerthoud et al.l ( 120071) . who fit 
the CO P' overtone emission of the Be star 51 Oph with a circum- 
stellar disc model. The line widths are also approximately ten times 
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Table 3. The parameters of the best-fitting models. The mass is determined from the luminosity presented in Table[T] The uncertainties in i, A'COj,, and Av are 
determined by holding all other parameters constant and determining the change required to increase the reduced ^^ by { q 
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0.59 


5.08 


2.2+5-'xl02i 


2.79 
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0.45 
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1.7+^^x1022 
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0.62 



i: No uncertainties are presented as these inner radii are the inner radii of the region where CO emission is possible, i.e. T< 5000K. Therefore, decreasing the 

inner radii did not affect the quality of the fits. 

*: No uncertainties are presented as these outer radii are the outer radii of the region where CO emission is possible, i.e. T> lOOOK. Therefore, increasing the 

outer radii did not affect the quality of the fits. 

t: ^2 flat down to 10"*cm-2, the minimum surface number density considered. Therefore, no error is presented here. 



the thermal broadening due to motion of the CO molecules. There- 
fore, the dominant contribution must be due to turbulence. How- 
ever, we do not consider shear broadening in the model. If this were 
taken into account, the turbulent velocity may well be less than the 
best-fitting values presented in Table[3] The properties of accretion 
discs around MYSOs ha s only recently begun to be examined in 
detail jVaidva et alj2009f) , and incorporating vertical disc structure 
to the disc model (which is required to assess sh ear broadening 
iHome & Marsi]|l986l : JHummel & Vranckerj|2000l) is beyond the 
scope of this paper. Therefore, we neglect shear velocity. However, 
it is important to note that neglecting vertical structure may arti- 
ficially limit the extent of the CO emitting region. In a disc with 
a vertical extent the upper regions of the disc may be hot enough 
to excite CO emission at radii where the mid-plane temperature 
has dropped below the required ~ 2000 K. Therefore, incorporat- 
ing vertical structure in the model may allow CO emission from 
a wider range of radii than the tliin disc model. This would have 
the effect of increasing the resulting spectro-astrometric signatures, 
and thus the predictions presented here may be lower limits. 



4.2 Comparisons with previous work 

Little is known about the circumstellar environments of the RMS 
objects. In particular, there are no previous high resolution stud- 
ies with which to compare the best-fitting model parameters. How- 
ever, the two non-RMS objects in the sample, M8 E and IRAS 
08576-4334 have been both been studied previously. Here, we as- 
sess whether the best-fitting models for these objects are consistent 
with other observations. 

iBik & Thil ( l2004h determine the inclination of IRAS 
08576-4334 to be 2T\l based on fitting the CO P' overtone band- 
head. The best-fitting value of ~ 18° determined here is thus con- 
sistent with the previous value. This might be expected as we apply 
the same methodology, none the less, this provides an important 
check th at the results are co nsistent with previous work. Turning 
to M8E, ISimon et alj l ll985l) postulate that this object possesses 
an edge o n disc; contrary t o our best-fitting inclination of ~ 16°. 
However, iLinz et al.l 1 120091) find an inclination < 30° is required 
to fit the SED of M8E. These auth ors suggest that the interpreta- 
tion of the lunar occultation data of ISimon et alj jl985r) is compli- 
cated by the presence of scattered l ight a nd outflow cavities. Given 
that the conclusion of iLinz et al.l 1 120091) is based on the SED of 
M8 E from the visible to mid infrared, in addition to the grid of 



models of iRobitaille et alj 1 2007D , we suggest the inclination de- 
termined bv lLinz et al.l ( 120091) is currently the best estimate, which 
is in agreement with the best-fitting value. Therefore, we conclude 
that the best-fitting models are consistent with results in the liter- 
ature, although only a few are available, indicating that the best- 
fitting model parameters are representative of the circumstellar en- 
vironments of the sample. 



4.3 G287.3716-h00.6444 

In all but one case, the disc model not only provided a good fit to 
the data (as measured with x^), but was also consistent with the 
observed flux densities, and where they existed, previous observa- 
tions. Therefore, it would appear that small scale discs are present 
around the majority of the sample. However, it was not possible to 
fit the CO P' overtone bandhead of G287.3716+00.6444 with the 
circumstellar disc model used to fit the other profiles. As shown in 
Fig.[2]the disc model could not fit observed the bandhead shoulder, 
nor the slope red-wards of the peak. Furthermore, the bandhead 
could not be fit with emission from an isothermal, non-rotating 
body of CO. Therefore, we are led to consider alternative scenarios 
to explain the emission and its profile. 

Besides hot, dense discs there are other viable sources of CO 
P' overto ne bandhead em i ssion. One such scenario is a dense, neu- 
tral wind. IChandler et al.l ( Il995l) were able to fit the CO P' over- 
tone bandheads of several YSOs with models of neutral winds, 
but note that th e mass loss rate s requi red were relatively high, up 
to 10'^ Mnvr~' . IChandler et akl ( Il995l) consider it unlikely CO P' 
overtone emission originates in a wind as such mass loss rates are 
much higher than observed for solar mass YSOs. However, the 
winds of MYSOs may well lead to mass loss rates of 10"^ Moyr"' 
JDrew et al.ll 19930 . Consequently, CO P' overtone emission from 
a dense wi nd should perhaps b e re-considered in this case. Al- 
ternatively, [Scoville_et_alJ ([l983]) propose that the CO P' overtone 
bandhead emission of the Becklin-Neugebauer object is created in 
shocks (as the observed velocity dispersion and estimated emitting 
area are both small). 

We note that iBlum et al.l ( |2004|) also report that the CO P' 
overtone bandhead of one of their sample (of four YSOs) was dif- 
ficult to fit with a model of CO emission fro m a circumstellar disc. 
Following the example of lKraus et al.l 1 120001) , they suggest that this 
may be due to the circumstellar disc exhibiting an outer bulge, 
which shields the inner regions of the disc. The effect of this would 
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be to limit tiie visible CO emitting region to low velocity regions, 
resulting in a narrow profile. In this case, however, it is not the 
width of the profile we cannot fit, but rather the slope of the profile 
red-wards of the bandhead peak. Most formation scenarios, such as 
winds and discs, result in excess blue-shifted emission, therefore 
this bandhead profile is difficult to explain. 

It is conceivable that the emission from several, discrete re- 
gions in a shock, will have different excitation temperatures, and 
when superimposed, will result in a different slope to the bandhead 
than the disc model. If such shocks exist, the shocked region must 
be located close to the central star, within a few au, as we do not 
see a positional excursion in the spectroastrometric data. However, 
it is difficult to envisage a scenario in which the shock emission 
is predominately red-shifted. As an alternative it may be that the 
emission does originate in a disc, but the receding part of the disc 
is significantly brighter than the approac hing side, similar to th e 
V/R variations exhibited by Be stars (e.g. iHanuschik et alJI 19950 . 
Regardless, it is apparent that while the majority of MYSO CO 
bandheads are well fit by models of circumstellar discs, the circum- 
stellar environments of MYSOs are not yet completely understood. 
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5 CONCLUSIONS 

In this paper we present high resolution near infrared spectroas- 
trometry around 2.3yum of a sample of young stellar objects, most 
of which are drawn from the RMS survey. The RMS constitutes 
the largest catalogue of massive young stellar objects to date, and 
thus provides a unique sample to study massive star formation. We 
model the CO P' overtone bandhead emission detected with a sim- 
ple model of a circumstellar disc in Keplerian rotation. In addition, 
the sub milli-arcsecond precision spectroastrometric data are com- 
pared to the spatial signatures of the best-fitting models. 

The observed bandhead profiles are, on the whole, well fit by 
models of the emission originating in circumstellar discs. This con- 
curs with the generally accepted view that this emission arises in 
small scale accretion discs. The observed bandhead of one object 
cannot be fit with a circumstellar disc model. This could be due to 
the emission emanating from an asymmetric disc or shock. This 
highlights that the circumstellar environments of massive young 
stellar objects are still not completely understood. 

No spatial signatures of discs are revealed in the spectroastro- 
metric data, which have a precision of approximately 0.4 mas. This 
is entirely consistent with the sizes of the best-fitting models. Due 
to the small sizes of the discs, the large distances to them and the 
decrease in brightness with radius the predicted signatures are of 
the order of O.Olmas. This is well below the current detection limit. 
However, the predicted signatures could perhaps be revealed by dif- 
ferential phase observations with AMBER. Currently, MYSOs are 
below the sensitivity limit of AMBER, but PRIMA should allow us 
to observe such objects and probe the spatial distribution of the CO 
emission at an even higher precision. 

To summarise, in general the model of emission from a cir- 
cumstellar disc provides a good fit to the observed bandheads, and 
is consistent with the observed flux densities. This indicates that the 
majority of MYSOs with CO P' overtone emission possess small- 
scale, circumstellar discs. In turn, this provides further evidence 
that massi ve stars form via disc ac cretion, as suggested by the sim- 
ulations of lKrumholz et al.l ( 120090 . 
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